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We investigate the effect of ferroelectric polarization direction on the geometric properties of Pd 
deposited on the positive and negative surfaces of LiNbOs (0001). We predict preferred geometries 
and diffusion properties of small Pd clusters using density functional theory, and use these calcu- 
lations as the basis for kinetic Monte Carlo simulations of Pd deposition on a larger scale. Our 
results show that on the positive surface, Pd atoms favor a clustered configuration, while on the 
negative surface, Pd atoms are adsorbed in a more dispersed pattern due to suppression of diffusion 
and agglomeration. This suggests that the effect of LiNbOs polarization direction on the catalytic 
activity of Pd [J. Phys. Chem. 88, 1148 (1984)] is due, at least in part, to differences in adsorp- 
tion geometry. Further investigations using these methods can aid the search for catalysts whose 
activities switch reversibly with the polarization of their ferroelectric substrates. 
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The spontaneous polarization of ferroelectric materi- 
als has enabled their use in technological devices ranging 
from SONAR to random access memory. While most cur- 
rent applications of ferroelectrics result primarily from 
their bulk properties, polarization can also impact the 
surface properties of these materials, including surface 
stoichiometry [1-3], geometry [2, 3], and electronic struc- 
ture [1, 3-5]. These polarization dependent differences in 
intrinsic surface properties also affect their interactions 
with adsorbates, [6-13] as evidenced by differences in ad- 
sorption energies [6-10] or the rates of surface catalyzed 
reactions [11-13]. Polarization orientation can, in turn, 
affect the chemical properties of the adsorbates them- 
selves. For example, Inoue et al. [11] showed that the 
activation barrier for CO oxidation by Pd adsorbed on 
a LiNbOs surface changes by 30 kJ/mol, depending on 
polarization orientation. The notion that a ferroelectric 
substrate's polarization can affect the activity of a sup- 
ported catalyst suggests the intriguing possibility that 
the activity of a catalyst could be modulated reversibly 
by switching the polarization of a ferroelectric substrate. 

Despite ample evidence showing that metals on oppo- 
sitely poled ferroelectric surfaces have different catalytic 
properties, the mechanism underlying these differences is 
not well understood. The most prevalent explanation for 
this phenomenon has been that the difference in charge 
between the two surfaces alters the electronic proper- 
ties of the catalyst [4, 11-13]. However, this explana- 
tion may be incomplete, as metal adsorption geometries 
significantly impact their catalytic properties [14]. Since 
ferroelectric surfaces may have different geometries de- 
pending on the sign of their polarization [1-3], it is plau- 
sible that the geometries of metals adsorbed onto them 
also differ [10, 15]. Two recent studies of Pd adsorp- 
tion on oppositely poled LiNbOs (0001) surfaces present 
conflicting conclusions regarding the possible effect of po- 



larization on metal adsorption geometry. Yun et al. [15] 
showed that large Pd clusters form on both the positive 
(c + ) and negative (c~) surfaces of LiNbOs, suggesting 
that polarization has little impact on metal adsorption 
geometry. In contrast, Zhao et al [10] observed large Pd 
clusters only on the c + surface, and a more planar ge- 
ometry on the c~ surface, suggesting that polarization 
strongly affects metal adsorption geometry. In addition, 
only Zhao et al observed a difference in CO tempera- 
ture programmed deposition between the two surfaces, 
suggesting that when polarization affects the activity of 
adsorbed catalysts, it does so, at least in part, by altering 
their geometries. 

In view of the ambiguity of experimental investigations 
of the relationship between polarization and metal ad- 
sorption geometry, we address this question on a micro- 
scopic level using theoretical methods. Here we inves- 
tigate the energetics and kinetics of Pd deposition on 
LiNbOs (0001) surfaces using a combination of density 
functional theory (DFT) calculations and kinetic Monte 
Carlo (KMC) simulations. We first calculate the adsorp- 
tion geometries of clusters on the c + and c~ surfaces. We 
then model the range of possible diffusion and agglomer- 
ation processes of these clusters using the nudged elastic 
band (NEB) method [16]. Finally, we use the activation 
barriers of these processes as inputs for a KMC simula- 
tion of the deposition of Pd on LiNbOs on a larger scale. 
To our knowledge, this is the first theoretical study of 
metal adsorption kinetics on a ferroelectric surface. 

In all DFT calculations, we studied five trilayer thick 
LiNbOs slabs with y/3 x y/3 surface supercells (Supple- 
mentary Fig. SI). We passivated surface charges with 
one Li atom per primitive supercell on the c + surface, 
and one O and one Li atom on c~, in accordance with 
the findings of Levchenko and Rappe that this is the 
most thermodynamically stable surface composition for 
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FIG. 1: (Color online) (a) Left Side: Monomeric potential 
energy surface and minimum energy geometries of Pdi, Pd2, 
Pd3 and Pd4 on c + . Contour spacing is 0.1 eV. Right side: 
Minimum energy geometries of (b) Pdi, (c) Pd2, and (d) Pds 
on c~ . Triangles are original positions of Li' . The geometries 
for Pd3 and Pd4 are similar to those of atoms 1-2-3 and 1-2- 
3-4 in (d). 



LiNb03 [2]. Hereafter, we denote these passivation atoms 
as Li' and O'. To remove spurious interactions between 
slabs in different supercells, we separated the slabs with 
more than 13 A of vacuum and employed a dipole cor- 
rection [17]. DFT total energy calculations [18] were 
performed using the generalized gradient approximation 
(GGA-PBE) [19] and implemented using the Quantum 
Espresso package [20]. Atoms were represented using 
norm-conserving nonlocal pseudopotentials [21, 22] gen- 
erated using the OPIUM code [23]. We allowed relax- 
ation of Pd atoms and the first two layers of the LiNbOa 
surface, but kept the remaining portion of the slab fixed. 

Our DFT results for Pd adsorption show differences in 
preferred binding geometries and diffusion barriers be- 
tween the c + and c~ surfaces. On the c + surface, Pd ad- 
sorption minimally changes the geometry of the LiNbOa 
surface itself. As a result, the monomeric potential en- 
ergy surface (PES, Fig. 1(a)) essentially determines the 
geometries of multiple Pd atoms adsorbed on this sur- 
face. In particular, all Pd atoms that we predict to bind 
directly to the c + surface prefer sites close to monomeric 
potential energy minima. However, because the distances 
between minima do not match well with optimal Pd-Pd 
bond lengths, binding of large numbers of Pd atoms both 
to the c + surface and to each other is unfavorable. For ex- 
ample, adsorption geometries of three and four Pd atoms 
both include one atom that bonds only to other Pd atoms 
and does not interact directly with the surface (Fig. 1(a)). 
Further, we found no metastable planar structure of four 
atoms, suggesting that the Pd agglomeration barrier on 
the c + surface is negligible. 

Diffusion and agglomeration processes on the c + sur- 
face are also impacted by the monomeric PES. In the 
monomeric case, this relationship is direct, as the two 
unique paths for monomer hopping (Fig. 2(a)) are the 
only paths containing saddle points on the PES. Be- 



cause the diffusion processes of clusters of up to four 
atoms on the c + surface are dominated by the move- 
ment of a single atom, they have diffusion barriers similar 
to monomer hopping, a fact that can also be explained 
by the monomeric PES. For example, in dimer walking, 
a process in which one atom in a dimer steps between 
two potential minima and stays bonded to the other Pd, 
paths 2 and 3 have activation barriers similar to monomer 
hopping paths 1 and 2, respectively (Fig. 2(b)). Sim- 
ilarly, though dimer sliding involves movement of both 
Pd atoms, one stays adjacent to the same Nb atom and 
thus, its movement contributes little to the activation 
barrier (Fig. 2(c)). Finally, because tetramer rolling re- 
quires movement of only one Pd atom out of a monomeric 
potential well, it too has an activation barrier similar to 
that of monomer hopping. The combination of low Pd 
diffusion barriers and a negligible agglomeration barrier 
suggests that formation of large Pd clusters is favorable 
on the c + surface. 

In contrast to its behavior on the c + surface, Pd ad- 
sorption on the c~ surface substantially alters the ge- 
ometry of the surface itself. This is largely due to in- 
teractions with O' atoms, which along with Li' atoms 
terminate the c~ surface for charge passivation [2]. The 
formation of Pd-O' bonds makes adsorption onto the c~ 
surface much more favorable than on c + (Supplementary 
Table SI). However this also leads to complex adsorption 
geometries, as the O' atom is free to tilt its bond to Nb 
substantially in order to accommodate bonding to Pd. 
Despite this geometric complexity, we find that the Pd 
adsorption geometries on the c~ surface can be under- 
stood in the context of maximizing the number of Pd-Pd 
and Pd-O' bonds formed. For example, this explains why 
the 2D planar (not shown) and 3D configurations of clus- 
ters of four or five Pd atoms on the c~ surface, which have 
identical numbers of Pd-Pd/Pd-O' bonds, have binding 
energies within 0.1 eV of one another. 

The inherent assumption in this analysis, that Pd-O' 
and Pd-Pd bonds have similar strengths, is justified by 
two facts. First, the average energies of these bonds, de- 
fined as adsorption energies divided by the number of 
Pd-Pd and Pd-O' bonds, in clusters of 1-5 Pd atoms are 
uniform (1.03 ± 0.05 eV). Second, these average bond 
energy values are similar to bond strengths in free clus- 
ters of 2-5 Pd atoms (1.06 ± 0.08 eV). From here on, we 
denote Pd-Pd and Pd-O' bonds as Pd-X. 

Just as Pd-X bonds are the primary contributors to 
the favorable adsorption energy of Pd clusters on the c~ 
surface, breaking these bonds, particularly Pd-O', is the 
primary barrier to cluster diffusion and agglomeration. 
Based on the argument above that Pd-X bonds have sim- 
ilar energies, we would expect the activation barriers for 
processes that require one Pd-O' bond breaking to be 
approximately 1 eV. However, actual activation barriers 
are lower, because other atoms can move to more favor- 
able positions in the process of transition state formation. 
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FIG. 2: (Color online) Schematic drawings of kinetic events. 
Top: Diffusion events on c + . (a) Pdi hopping, (b) Pd2 walk- 
ing, (c) Pd2 sliding and (d) Pd4 rolling. Bottom: Diffusion 
events on c~ . (e) Pdi hopping, (f) Pd4 rolling, (g) Pd4 ag- 
glomeration, and (h) agglomeration at large cluster (path 1) 
and in-plane motion (path 2). Large black circles and open 
blue circles denote initial and final positions for diffusion 
events, after the events. Blue dashed lines in (e) indicate 
that the monomer hops between O' bridge sites. 



This is especially true when new Pd-X bonds are formed 
before a Pd-X bond is fully broken. Consistent with our 
understanding of the effect of Pd-X bond breaking on the 
activation barriers of diffusion processes, we find that in- 
plane movement at the boundary of a cluster has a high 
activation barrier (path 2 in Fig. 2(h)), because it re- 
quires breaking of two Pd-X bonds. 

Our DFT calculations show that diffusion barriers are 
substantially lower on the c + surface (< 0.4 eV) than on 
the c~ surface (> 0.8 eV). These activation energies, E a , 
can be compared in a physically meaningful way when 
converted to expected event-event time intervals, r, us- 
ing Arrhenius kinetics [24], - = ve ~ E a/k B T ^ Assum- 
ing an attempt frequency (y — 10 12 sec -1 [25]) the time 
scales of diffusion events on the c + and c~ surfaces are 
on the order of microseconds and minutes, respectively. 
We can then infer that at a deposition rate of ~ 0.01- 
0.1 ML/s [15], on average, each new monomer deposited 
on the c + surface will aggregate to an existing cluster 
before the next atom is deposited. In contrast, we would 
expect many Pd atoms to be deposited in the vicinity 
of a given Pd atom on the c~ surface between diffusion 
events of that atom. We thus infer that Pd atoms will 
agglomerate into much larger clusters on the c + surface 
than on the c~ surface of LiNbOa. 

We conducted kinetic Monte Carlo (KMC) simulations 
in order to characterize this inferred difference in adsorp- 
tion geometries. Onto a 10 nm x 10 nm surface, we ran- 
domly deposited Pd atoms one by one at a selected de- 
position rate. We then allowed Pd atoms and clusters to 
attempt a series of diffusion and agglomeration processes, 
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TABLE I: Diffusion activation barriers (in eV) of processes 
drawn in Fig. 2 and described in the text. Barriers of re- 
verse processes are written below forward processes (c + ) or 
in parentheses (c~). 



with probabilities based on our calculated activation bar- 
riers, at a constant attempt frequency of 10 12 sec -1 [26] 
for all events (Full KMC description in Supplement). 

Our KMC simulations confirm that Pd forms larger 
clusters on the c + surface than on the c~ surface. Corre- 
spondingly, we find that Pd covers a much smaller area 
of the c + surface than of c~ . We also find that Pd area 
is insensitive to deposition rate at room temperature, as 
all diffusion and agglomeration processes are highly ac- 
tivated (for c + ) or suppressed (for c~) at this tempera- 
ture. On the c~ surface, the area covered by Pd is only 
affected when agglomeration plays a role, which occurs 
when coverage is above 0.5 ML (Fig. 3(a)). On the c + 
surface, Pd coverage area increases as temperature is low- 
ered (Fig. 3(b)). 

The relatively large diffusion barriers on the c~ sur- 
face (~ 0.8 eV) were not sufficient to prevent agglomera- 
tion completely, especially when our simulations were ex- 
tended until slightly after deposition was complete. How- 
ever, we infer that the clusters on the c~ surface remain 
much smaller than those on the c + surface over a very 
long time scale (Supplementary Fig. S3). 

Overall, our results support the conclusion of Zhao et 
al. that Pd cluster sizes are much larger on the c + sur- 
face than on c~ [10]. However, Yun et al. also produced 
a thorough data set supporting their conclusion that Pd 
forms large clusters on both surfaces of LiNbOa [15]. One 
possible explanation is that Yun and Zhao studied Pd 
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FIG. 3: (Color online) Proportion of surface area covered by 
Pd for simulations at a range of (a) deposition rates at 300 K 
and (b) temperatures with a deposition rate of 0.025 ML/s. 
KMC snapshots of 1.0 ML on (c) c + and (d) c~ at 300 K with 
a 0.025 ML/s deposition rate. 



adsorption onto LiNb03 surfaces with substantially dif- 
ferent compositions. This could be due to the presence of 
different impurities, such as hydroxyl groups [27], which 
are known to exist on some LiNbOs surfaces, or oxy- 
gen lattice vacancies [3]. Our calculations predict that 
Pd binding is weaker by 0.6 eV to a c~ surface termi- 
nated with OH instead of OLi, suggesting that agglom- 
eration would be more favorable with some OH impuri- 
ties present. Thus, it is plausible that if there were OH 
impurities on the c~ surface studied by Yun et al, they 
may have influenced the formation of larger Pd clusters. 

In conclusion, we find that the different surface geome- 
tries of oppositely poled LiNbOs (0001) surfaces lead to 
substantially different Pd adsorption geometries. On the 
c~ surface, strong Pd bonding to O', the extra surface 
oxygen present for charge passivation, leads to larger dif- 
fusion and agglomeration barriers. This lead in turn to 
less clustering and a more planar geometry overall on 
this surface than on c + . In contrast to the prevalent 
view that differences in catalyst charging are responsible 
for differences in the catalytic properties of metals ad- 
sorbed on polar ferroelectric surfaces, we conclude that 
the difference in adsorption geometry predicted here is 
sufficient to explain much of the difference in catalytic 
activity that has been observed for Pd deposited on op- 
positely poled LiNbOs surfaces [11]. Because formation 
of large Pd clusters is thermodynamically favorable on 
both c + and c~, it is unlikely that the catalytic activ- 
ity of Pd could be switched reversibly by switching the 
polarization of the LiNbOs substrate. However, the com- 
bination of theoretical methods used here is well suited 
for the study of other metal/ferroelectric surface com- 
binations in search of systems with switchable catalytic 
activity, and the present example serves as a paradigm 
of polarization controlling catalytic cluster geometry and 
reactivity. 
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ATOMIC MODEL OF LITHIUM NIOBATE 
SUBSTRATES 



The ferroelectric phase of lithium niobate (LiNbOs) 
has a bulk structure consisting of layers in a Li-03-Nb 
pattern [2]. Here we show LiNbOs with its thermody- 
namically preferred surface terminations, as predicted by 
Levchenko and Rappe (Fig. SI). Li atoms are added to 
the positive (c + ) surface, and both Li and O atoms are 
added to the negative surface (c~ ) of LiNbOs fc> r charge 
passivation. Passivation atoms are denoted Li 7 and O'. 
These surface terminations make the overall stoichiome- 
try of five trilayer slab used in this study, from c + (left) 
to c~ (right), Li / -(Li-0 3 -Nb) 5 -Li / / . 
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TABLE SI: Per- atom adsorption energies of Pd on LiNbOs 
(LNO), relative to infinitely separated Pd atoms (E^ ds ) and 
free Pd clusters (E^ ds ). We define E^ ds — (^pcI-lno-^lno- 
nE Pd )/n and E c ads = (^pd-LNO - -Elno - -Epd n )/n, where n is 
the number of Pd atoms. On the c~ surface, we report average 
Pd-X bond energies, E™~*= riE^ ds /N hond , where iVbond is 
the number of Pd-X bonds. We also report average Pd-Pd 
bond energies, ^bond Pd > m ^ ree clusters for comparison. 
All energies are given in eV. 



SUPPLEMENTARY DFT METHODS 

In all density functional theory (DFT) calculations, 
we used 3x3x1 /c-point grids [3], plane wave energy 
cutoffs of 50 Ry [4], and force tolerances of 0.01 eV/A. 
When necessary, we included spin as a degree of free- 
dom in our calculations. We checked the convergence of 
our monomer and dimer adsorption energies calculated 
in y/3 x y/3 supercells on five layer thick slabs by using 
either y/7 x \fl supercells or seven trilayer thick slabs. 
The resulting energies agreed within 0.1 eV/Pd. 
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FIG. SI: Geometry of five-trilayer LiNbOs slab: (a) side view, 
(b,c) top view of c + , and (d,e) top view of c~ . Slabs are 
represented using ball and stick (a) and space filling (b,d) 
models, or the simplified model used in the main text (c,e). 
Blue arrows in (a) denote ions (Li' and O') introduced for 
surface charge passivation. Black lines in (b) and (d) indicate 
the boundaries of the LiNb03 primitive surface unit cell, (b) 
and (d) are redrawn from Ref. [1]. 



ADSORPTION ENERGIES 

In our analysis of the energetics of Pd adsorption on 
LiNbOs, we consider two quantities for per-atom adsorp- 
tion energy, E% ds and E° ds , defined as the difference be- 
tween the energy of Pd-LNO and the sum of the energies 
of the LiNbOs slab and either infinitely separated (E% d3 ) 
or clustered (E^ ds ) Pd atoms divided by the number of 
Pd atoms (Table SI). We note that on the c + surface, 
E° ds becomes very small as the number of adsorbed Pd 
atoms increases, because of a mismatch between ideal 
Pd-Pd bond lengths and distances between Pd binding 
sites, and because most of the cluster not in contact with 
LNO. In contrast, E° ds remains relatively high even for 
pentamers on c~, largely due to Pd-O 7 bonding. In the 
main text, we explained many of our conclusions about 
energetics on the c~ surface in terms of the total num- 
ber of Pd-Pd and Pd-O 7 bonds (hereafter denoted Pd- 
X bonds). This analysis was justified by the fact that 
average Pd-X bond energies (^bond X )> defined as nE® ds 
divided by the number of Pd-X bonds (TVbond), were sim- 
ilar for all cluster sizes studied and were also similar to 
average bond energies in free Pd clusters (£^ond)- ^ e 
show explicit values for ^b d nd Pd and in Table SL 

Since Pd-O 7 bonding makes the key contribution to the 
difference in adsorption behavior between the c + and c~ 
surfaces, we wanted to be sure that DFT-GGA does not 
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FIG. S2: An FCC lattice is not sufficient to allow all possible 
agglomeration events to occur. Of the two agglomeration pro- 
cesses pictured, one is allowed by the FCC lattice, while the 
other, denoted by a red 'X', would also require HCP lattice 
points on the second layer of the simulation lattice. We thus 
used an A-(ABC)-(ABC)- • • lattice as described in the text. 
Open circles denote Pd atoms and 'X's denote points on the 
second layer of the simulation lattice. 

misestimate the favor ability of Pd-0 interactions. To do 
this, we considered the case of bulk PdO. We compared 
our DFT-GGA energy of formation of PdO from individ- 
ual Pd and O atoms, defined as E^ om + E^ m - E Pd0 , 
to the sum of the experimental values for PdO formation 
energy [5], Pd bulk cohesive energy [6], and half of the 
O2 bonding energy [7]. Our DFT-calculated value was 
only 2.1% greater than the experimental value, suggest- 
ing that DFT predicts Pd-0 interactions well and that 
the strong Pd-0 interactions seen on the c~ surface are 
valid. As a note, we did not use the DFT PdO formation 
energy as a metric for comparison, because DFT-GGA is 
known to predict the energy of O2 quite poorly [7]. 



KINETIC MONTE CARLO SIMULATION 
DETAILS 

Our kinetic Monte Carlo (KMC) simulations were per- 
formed in periodic 10 nm x 10 nm cells. We included 15 
layers of Pd binding sites above the surface, though no Pd 
atoms reached the top layer of the cell in any of our simu- 
lations. All Pd binding sites on the layer in contact with 
the LiNbOa surface were mapped onto their correspond- 
ing points on the triangular oxygen lattice. Beginning 
with the second layer above the surface, we used three 
sublattices to allow all possible Pd cluster configurations 
such that the lattice is substrate- A-( ABC)-(ABC)- •• . 
The use of this lattice, which is much denser than the 
face-centered cubic (FCC) lattice of bulk Pd, was neces- 
sary to describe some diffusion events that can occur for 
small numbers of Pd atoms, but is forbidden by the FCC 
lattice. For example, half of all possible agglomeration 
processes are forbidden using an FCC lattice, because 



lattice points are present that correspond to only half of 
the surface's hollow sites (Fig. S2). Additionally, hopping 
of Pd atoms between FCC and hexagonal close-packed 
(HCP) hollow sites on (111) cluster facets is forbidden 
by a pure FCC lattice. To compensate for the fact that 
some points in the lattice used in our simulation are closer 
together than a reasonable Pd-Pd bond distance, we im- 
posed two rules. First, the same site cannot be occupied 
on adjacent layers of the lattice. Additionally, within the 
same layer, adjacent sites on different sublattices cannot 
be occupied. 

We built our simulations primarily based on our DFT 
calculations of the activation barriers of diffusion and ag- 
glomeration of up to four atoms. On the c + surface, it 
was unnecessary to consider diffusion of larger clusters, 
because our deposition time interval was sufficiently long 
that we observed monomers to aggregate into clusters be- 
fore other monomers are deposited in the same vicinity. 
Diffusion of large clusters on the c~ surface could also 
be neglected, because they require the breaking of two 
or more Pd-X bonds and thus have activation barriers 
sufficient to prevent their occurrence on the time scale 
that we considered. As described below, in cases where 
we had not calculated activation barriers directly, such 
as monoatomic motion on cluster facets, we interpolated 
reasonable barriers. 

We permitted diffusion of Pd atoms that were directly 
in contact with the LiNbOa surface or on cluster facets. 
On the LiNbOa surface, we permitted the following diffu- 
sion processes: 1) nearest-neighbor monomer hopping, 2) 
agglomeration/deagglomeration (movement of Pd atom 
from layer 1 (2) onto layer 2 (1)), 3) dimer walking and 
sliding on c + , 4) trimer walking and sliding on c + , 5) 
trimer flipping on c~, 6) tetramer rolling and 7) dissoci- 
ation (some of these processes are pictured in Fig. 2 in 
the main text). On cluster facets, we permitted monomer 
hopping only. Though concerted movements of multiple 
atoms can occur on metal cluster facets, these generally 
affect cluster morphology, but not cluster size (number 
of atoms present). Since our primary interest is only in 
seeing differences in cluster size, this simplification of Pd 
behavior on cluster facets is acceptable. 

It is well known that the activation barriers for 
monomeric diffusion on the (111) surfaces of FCC metals 
are quite low relative to the bulk metal-metal bond en- 
ergy. Our nudged elastic band (NEB) calculations show 
that the activation barrier for monoatomic hopping from 
from FCC to HCP hollow sites on (111) surface of pal- 
ladium is 0.11 eV, despite the fact that the Pd-Pd bulk 
bond energy is nearly six times larger (0.59 eV). This 
is explained by the fact that during monomer hopping 
between nearest-neighbor sites, new bonds start to form 
before others are completely broken. In our KMC sim- 
ulations, we interpolated reasonable activation barriers 
for diffusion events on Pd cluster facets using insights 
from Trunshin and coworker's studies of nearest-neighbor 
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FIG. S3: Pd adsorption geometries predicted by kinetic Monte Carlo simulations on the c~ surface of LiNbOs immediately 
after deposition ends (a-c) and 2-15 minutes later (d-e). These simulations were conducted at different coverages (0.1 ML 
and 1.0 ML), different temperatures (300 and 400 K), and with the two different agglomeration barriers predicted by DFT 
(0.59 eV and 0.74 eV). All simulations were conducted at a deposition rate of 0.025 ML/s. Times waited after the completion 
of deposition were: (d) 15 minutes, (e) 10 minutes, (f) 2 minutes, and (g) 3 minutes. All simulations are conducted at 300 K, 
except for (f), which is at 400 K. 



monoatomic hopping on copper surfaces and steps [8]. 
By compensating for the difference in cohesive energy 
between bulk Pd and Cu, we estimated that the acti- 
vation barrier for nearest neighbor monomer hopping is 
0.2 eV per Pd-Pd bond broken during transition state 
formation. We also showed that the qualitative results 
of some of our simulations were insensitive to the use 
of a 0.3 eV barrier for this process. In next nearest 
neighbor monomer hopping, which is analogous to the ag- 
glomeration processes pictured in Fig. 2 but generalized 
to movements between any two levels in the simulation 
lattice (starting with level 2), Pd-Pd bonds are broken 
completely during transition state formation. Thus, we 
estimated the activation barrier for this process based on 
the cost of Pd-Pd bond breaking. The cost of breaking 
a bond in a Pd2o cluster cut from bulk Pd is 0.52 eV. 
Similarly the cost of dimer dissociation is 0.55 eV. Thus, 
for next-nearest neighbor hopping on Pd cluster facets, 
we used a barrier of 0.5 eV per Pd-Pd bond broken. 



AGGREGATION AND AGGLOMERATION ON 
c~ AFTER DEPOSITION 

The relatively large diffusion barriers on the c~ surface 
(~0.8 eV) were not sufficiently high to prevent agglom- 
eration completely, especially when our simulations were 
extended until slightly after deposition was complete. We 
have conducted KMC simulations during deposition and 



after annealing for 2-15 minutes for both low (0.1 ML) 
and high (1.0 ML) coverages and at different temper- 
atures (300 or 400 K). We also conducted our simula- 
tions with agglomeration barriers (0.59 and 0.74 eV) for 
the two different possible cluster agglomeration processes 
predicted by DFT on the c~ surface (See Fig. 2(g,h) and 
Table I in the main text). These data are reported in 
Fig. S3. 

At low coverage (0.1 ML), Pd atoms are highly dis- 
persed at the end of deposition (Fig. S3 (a)). However, 
within 15 minutes, all but two Pd clusters are immo- 
bile, suggesting that on a slightly longer time scale, all 
Pd atoms will join immobile clusters (Fig. S3(d)). How- 
ever, these clusters are quite small compared to the sin- 
gle cluster that forms on the c + surface under the same 
conditions. Two factors contribute to the formation of 
smaller clusters on the c~ surface than on c + : lower max- 
imum mobile cluster size, and the higher number den- 
sity of clusters on the surface during deposition. On the 
c + surface, the mismatch between the distances between 
monomeric Pd binding sites and Pd-Pd equilibrium bond 
distances makes the activation barrier for motion of rel- 
atively large clusters rather low. For example, our pre- 
dicted barrier for tetramer rolling was similar to that for 
monomer hopping, largely because this process requires 
movement of only one Pd atom out of a monomeric po- 
tential well. Other work suggests that this condition may 
be met by Pd clusters as large as 13 atoms [9], and thus, 
that clusters of this size may be mobile on the c + sur- 
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face. In contrast, the diffusion of clusters larger than a 
tetramer on the c~ surface requires breakage of at least 
two Pd-X bonds, making these diffusion events extremely 
unlikely on the time scale of deposition. A second fac- 
tor contributing to the formation of large clusters on the 
c + surface is the fact that, because barriers for diffusion 
are so low, many diffusion events can occur during the 
average waiting time between addition of new atoms to 
the simulation cell on this surface. This means that only 
a few atoms at a time will be separate from the large 
immobile cluster, greatly reducing the probability that a 
second immobile cluster will form. The opposite is true 
on the c~ surface, where larger diffusion barriers mean 
that on average multiple Pd atoms are deposited inside 
the simulation cell between diffusion events of a single 
atom or cluster. As a result, many individual clusters 
are present in any given area of the simulation cell, pro- 
moting the formation of many small immobile clusters. 

At high coverage (1.0 ML), Pd atoms continue to re- 
arrange after deposition is complete until all atoms make 
at least three Pd-X bonds (Fig. S3(b,e)). Once this oc- 
curs, agglomeration can only occur on an extremely long 
time scale, because it requires breakage of multiple Pd-X 
bonds. Thus, the Pd coverage area of these clusters will 
remain relatively constant. This was true even when our 
surface was annealed at a higher temperature (400 K, 
Fig. S3(f)), consistent with the finding of Zhao and col- 
leagues that the Pd geometry on the c~ surface is stable 
up to 425 K [10]. 

Even when we used the lower of our two agglomera- 
tion barriers predicted from DFT on the c~ surface (Fig. 
2(g) in the main text), multiple distinct clusters formed 
at 1 ML coverage (Fig. S3(c,g)). At low coverage we 
observed results similar to those for the higher agglom- 



eration barrier (Fig. S3(d)), but with taller individual 
clusters. Overall, our simulations with this lower agglom- 
eration barrier on the c~ surface predicted Pd coverage 
areas that were substantially higher than those on the c + 
surface, but lower than the simulations on the c~ surface 
with the higher agglomeration barrier. 

In sum, our DFT and KMC calculations predict that 
though Pd forms slightly three dimensional clusters on 
the c~ surface of LiNbOs, its overall adsorption pattern 
on this surface is much more dispersed than that on the 
c + surface. 



corresponding author: rappe@sas.upenn.edu 
[1] S. V. Levchenko and A. M. Rappe, Phys. Rev. Lett. 100, 
256101 (2008). 

[2] A. Saito, H. Matsumoto, S. Ohnisi, M. Akai-Kasaya, 
Y. Kuwahara, and M. Aono, Jap. J. Appl. Phys 43, 2057 
(2004). 

[3] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 
(1976). 

[4] J. Ihm, A. Zunger, and M. L. Cohen, J. Phys. C 12, 4409 
(1979). 

[5] J. Nell and H. S. O'Neill, Geochimica et Cosmochimica 
Acta 60, 2487 (1996). 

[6] C. Kittel, Introduction to Solid State Physics (John Wi- 
ley & Sons, Inc., 1996), seventh ed. 

[7] F. Furche, Phys. Rev. B 64, 195120 (2001). 

[8] O. Trushin, A. Karim, A. Kara, and T. S. Rahman, Phys. 
Rev. B 72, 115401 (2005). 

[9] W. Fan and X. Gong, Applied Surface Science 219, 117 
(2003). 

[10] M. S. H. Zhao, D. A. Bonnell, and J. M. Vohs, J. Vac. 
Sci. Technol. 27, 1337 (2009). 



